theory BTO
imports Main

begin

(*

This file has two basic parts: a mathematical definition of matrix
operations and a BTO definition of the same matrix operations.  The
mathematical definition treats matrices as functions of type

(nat x nat) => nat

and composes them as per a typical mathematical textbook.  The BTO
version creates matrices as linked lists of linked lists, and each
operation builds a new matrix.  We have a well_formed predicate that
requires the linked lists to be of the same length (i.e. no ragged
matrices) and have sizes as expected.

Our compiler takes an AST and converts it into a series of Isabelle
function calls using the BTO matrix datatype.  We will prove this to be
equivalent to the pure mathematical definition of the same computations.

- Thomas Nelson and Jonathan Turner

*)

(*
-----
Types
-----
*)

types row = "nat list"
types matrix = "nat list list"

(* Note: Matrix contains the matrix and the size (which we use to check the matrix for well-formedness) *)
datatype expr = Matrix matrix nat nat | Add expr expr | Mult expr expr | Scale nat expr
datatype math_matrix = Mat "((nat \<times> nat) \<Rightarrow> nat)" nat nat 

(*
------
math matrix
------
*)

(* a convenience functions *)
primrec math_matrix_fun :: "math_matrix \<Rightarrow> ((nat \<times> nat) \<Rightarrow> nat)"
where
"math_matrix_fun (Mat f m n) = f"

(* mathematical addition *)
fun math_add :: "math_matrix \<Rightarrow> math_matrix \<Rightarrow> math_matrix"
where
"math_add (Mat A m n) (Mat B p q) = 
  (Mat (% i. (A i) + (B i)) m n)"
(* mathematical scalar multiplication *)
primrec math_scale :: "nat \<Rightarrow> math_matrix \<Rightarrow> math_matrix"
where
"math_scale alpha (Mat A m n) = (Mat (% i. alpha * (A i)) m n)"

(* mathematical matrix multiplication *)
fun math_mult :: "math_matrix \<Rightarrow> math_matrix \<Rightarrow> math_matrix"
where
"math_mult (Mat A m n) (Mat B p q) =
  (Mat (% (i,j). foldr plus (map (% k . A(i,k) * B(k,j)) [0..<n]) 0) m q)"

(*
---------
well-formedness
---------
*)

(* Well-formed row *)
inductive wf_row :: "row \<Rightarrow> nat \<Rightarrow> bool" where
wfr0[intro!]: "wf_row [] 0" |
wfr1[intro!]: "\<lbrakk>wf_row r1 n\<rbrakk> \<Longrightarrow> wf_row (e1#r1) (Suc n)"

(* Well-formed matrix *)
inductive wf_matrix :: "matrix \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> bool" where
wfmat0[intro!]: "wf_matrix [] 0 x" |
wfmat1[intro!]: "\<lbrakk>wf_row r1 row_len ; wf_matrix m1 n row_len \<rbrakk> \<Longrightarrow> 
	wf_matrix (r1#m1) (Suc n) row_len "

inductive wf_expression :: "expr \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> bool" where
wfexpr0[intro!]: "\<lbrakk> wf_matrix mat col_len row_len; col_len > 0; row_len > 0 \<rbrakk> \<Longrightarrow> wf_expression (Matrix mat col_len row_len) col_len row_len" |
wfexpr1[intro!]: "\<lbrakk> wf_expression e1 m n; wf_expression e2 m n\<rbrakk> \<Longrightarrow> wf_expression (Add e1 e2) m n" |
wfexpr2[intro!]: "\<lbrakk> wf_expression e1 m n; wf_expression e2 n p\<rbrakk> \<Longrightarrow> wf_expression (Mult e1 e2) m p" |
wfexpr3[intro!]: "\<lbrakk> wf_expression e m n \<rbrakk> \<Longrightarrow> wf_expression (Scale s e) m n"

inductive_cases wf_row_null: "wf_row r 0"
inductive_cases wf_row2[elim!]: "wf_row (x # xs) (Suc m)"
inductive_cases wf_row_succ: "wf_row r (Suc m)"
inductive_cases wf_row_empty: "wf_row [] m"

inductive_cases wf_matrix_null: "wf_matrix A 0 0"
inductive_cases wf_matrix_null2: "wf_matrix A 0 (Suc n)"
inductive_cases wf_matrix_null3: "wf_matrix A 0 n"
inductive_cases wf_matrix2[elim!]: "wf_matrix (x # xs) (Suc m) n"
inductive_cases wf_matrix_succ: "wf_matrix A (Suc m) n"
inductive_cases wf_matrix_empty: "wf_matrix [] n p"

(*
----------
BTO functions
----------
*)

(* BTO row addition *)
definition row_add :: "row \<Rightarrow> row \<Rightarrow> row" where
	"row_add x y = map (%(a,b) . a+b) (zip x y)"

(* BTO matrix add *)
definition matrix_add :: "matrix \<Rightarrow> matrix \<Rightarrow> matrix" where
	"matrix_add x y = map (%(a,b) . row_add a b) (zip x y)"

(* BTO matrix access (indexing operator) *)
definition matrix_access :: "matrix \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where
	"matrix_access m i j = ((m ! i) ! j)"

(* BTO row scale *)
definition row_scale :: "nat \<Rightarrow> row \<Rightarrow> row" where 
	"row_scale x = map (% a . a*x)"

(* BTO matrix scale *)
definition matrix_scale :: "nat \<Rightarrow> matrix \<Rightarrow> matrix" where
	"matrix_scale x = map (row_scale x)"

(* BTO vector matrix multiply *)
(* w' = u' B *)
fun vector_matrix_multiply :: "row \<Rightarrow> matrix \<Rightarrow> row"
where
	"vector_matrix_multiply [] [] = []" |
	"vector_matrix_multiply (e0 # u) (v0 # B) = 
		foldr (% (e,v) a . row_add (row_scale e v) a) (zip u B) (row_scale e0 v0)"

(* BTO matrix multiply *)
(* C = A B *)
definition matrix_multiply :: "matrix \<Rightarrow> matrix \<Rightarrow> matrix" where
	"matrix_multiply A B = 
		map (% u. vector_matrix_multiply u B) A"

(*
------
Computation & Evaluation
------
*)


(* BTO compute/evaluate *)
fun compute :: "expr \<Rightarrow> matrix" where
  "compute (Matrix m1 rows cols) = m1" |
  "compute (Add e1 e2) = matrix_add (compute e1) (compute e2)" |
  "compute (Scale s1 e1) = matrix_scale s1 (compute e1)" |
  "compute (Mult e1 e2) = matrix_multiply (compute e1) (compute e2)"

(* Mathematical evaluate *)
fun math_eval :: "expr \<Rightarrow> math_matrix " where
  "math_eval (Matrix m1 rows cols) = Mat (% (i,j). matrix_access m1 i j) rows cols" |
  "math_eval (Add e1 e2) = math_add (math_eval e1) (math_eval e2)" |
  "math_eval (Scale s1 e1) = math_scale s1 (math_eval e1)" |
  "math_eval (Mult e1 e2) =  math_mult (math_eval e1) (math_eval e2)" 

(* Show equivalence between the BTO computation and the mathematical model *)

(* 
a BTO matrix is equivelent to a math matrix iff they are the same at every index.

*)
inductive matrix_equivalence :: "math_matrix \<Rightarrow> matrix \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> bool" where
mat_eq[intro!]:"\<lbrakk> wf_matrix btomat m n ; (\<forall> i j . (i < m \<and> j < n) \<longrightarrow> (matrix_access btomat i j) = f (i,j)) \<rbrakk> \<Longrightarrow> matrix_equivalence (Mat f m n) btomat m n"

(* lifting a matrix into math_matrix preserves equivalence *)
lemma repr_equivalence: "wf_matrix A m n \<Longrightarrow> matrix_equivalence (Mat (%(i,j) . matrix_access A i j) m n) A m n"
apply (auto)
done

(* Utilities *)
lemma wf_row_implies_len: "wf_row r m \<Longrightarrow> length r = m"
apply (induct rule: wf_row.induct)
apply auto
done

lemma wf_matrix_implies_len: "wf_matrix A m n \<Longrightarrow> length A = m"
apply (induct rule: wf_matrix.induct)
apply auto
done

lemma row_len_in_matrix: "\<lbrakk> wf_matrix A m n; i < m \<rbrakk> \<Longrightarrow> wf_row (A ! i) n"
apply (induct arbitrary: i rule: wf_matrix.induct)
apply simp
proof -
	fix r1 row_len m1 n i
	assume 1: "wf_row r1 row_len" 
	and 2: "wf_matrix m1 n row_len" 
	and 3: "\<And>i. i < n \<Longrightarrow> wf_row (m1 ! i) row_len" 
	and 4: "i < Suc n"
	from 1 2 3 4 show "wf_row ((r1 # m1) ! i) row_len"
	proof (cases i)
		assume 5: "i = 0" 
		from 5 have 7: "(r1 # m1) ! i = r1" by simp
        from 7 1 show ?thesis by auto
	next
		fix nat
		assume 6: "i = Suc nat" 
		from 4 6 have "nat < n" by simp
		from this 3 have 8: "wf_row (m1 ! nat) row_len" by simp
		from 6 have "(m1 ! nat) = ((r1 # m1) ! i)" by simp
		from this 8 show ?thesis by simp
	qed
qed

(* Scale *)
lemma row_scale_wf: "wf_row r m \<Longrightarrow> wf_row (row_scale s r) m"
apply (induct arbitrary: s rule: wf_row.induct)
unfolding row_scale_def
apply auto
done

lemma matrix_scale_wf: "\<lbrakk> wf_matrix A m n\<rbrakk> \<Longrightarrow> wf_matrix (matrix_scale s A) m n"
apply (induct arbitrary: s rule: wf_matrix.induct)
unfolding matrix_scale_def
apply (auto simp: row_scale_wf)
done

lemma scale_row_access: "\<lbrakk> wf_row r m ; i < m \<rbrakk> \<Longrightarrow> (row_scale s r) ! i = s * (r ! i)"
unfolding row_scale_def
apply (drule wf_row_implies_len)
apply auto
done

lemma scale_access: "\<lbrakk> wf_matrix B m n; i < m; j < n \<rbrakk> \<Longrightarrow> matrix_access (matrix_scale s B) i j  = s * (matrix_access B i j)"
unfolding matrix_access_def
unfolding matrix_scale_def
apply (frule wf_matrix_implies_len)
apply (simp add: List.nth_map)
apply (drule row_len_in_matrix)
apply assumption
apply (erule scale_row_access)
apply assumption
done
	

lemma scale_equivalence: "matrix_equivalence mathM btoM m n \<Longrightarrow> matrix_equivalence (math_scale s mathM) (matrix_scale s btoM) m n"
apply (erule matrix_equivalence.cases)
apply auto
apply (force simp: matrix_scale_wf)
proof -
fix f m n i j
    assume 1:"wf_matrix btoM m n" 
    and 2:"\<forall>i j. i < m \<and> j < n \<longrightarrow> matrix_access btoM i j = f (i, j)" 
    and 3:"i < m" 
    and 4:"j < n"
	from 2 3 4 have 5: "matrix_access btoM i j = f (i, j)"  by simp
	from scale_access 5 1 3 4 show "matrix_access (matrix_scale s btoM) i j = s * f (i, j)" by simp
qed

(* add *)

lemma add_row_access: "\<lbrakk> wf_row r m ; wf_row s m; i < m \<rbrakk> \<Longrightarrow> (row_add r s) ! i = (r ! i) + (s ! i)"
unfolding row_add_def
apply (drule wf_row_implies_len)+
apply auto
done

lemma add_access: "\<lbrakk> wf_matrix A m n; wf_matrix B m n; i < m; j < n \<rbrakk> \<Longrightarrow> (matrix_access (matrix_add A B) i j) = (matrix_access A i j) + (matrix_access B i j)"
unfolding matrix_access_def
unfolding matrix_add_def
apply (frule_tac A="A" and m="m" and n="n" in wf_matrix_implies_len)
apply (frule_tac A="B" and m="m" and n="n" in wf_matrix_implies_len)
apply (simp add: List.nth_map)
apply (drule row_len_in_matrix)
apply assumption
apply (erule add_row_access)
apply (erule row_len_in_matrix)
apply auto
done

lemma wf_help: "\<lbrakk> wf_row r (Suc m) ; r = [] \<rbrakk> \<Longrightarrow> P"
apply (erule wf_row_succ)
apply blast
done

lemma row_add_wf: "\<lbrakk> wf_row r m; wf_row s m \<rbrakk> \<Longrightarrow> wf_row (row_add r s) m"
apply (induct m arbitrary: r s)
apply (erule wf_row_null)+
apply (simp add: row_add_def wfr0)
proof -
	fix m r s 
	assume 1: "\<And>r s. \<lbrakk>wf_row r m; wf_row s m\<rbrakk> \<Longrightarrow> wf_row (row_add r s) m" 
	and 2: "wf_row r (Suc m)" 
	and 3: "wf_row s (Suc m)"
	from 1 2 3 show "wf_row (row_add r s) (Suc m)" proof (cases r)
		assume 1: "wf_row r (Suc m)" and 2: "r = []" 
		from 1 2 show ?thesis by (rule wf_help)
	next
		fix a list
		assume 4: "r = a # list"
		from 1 2 3 4 show "wf_row (row_add r s) (Suc m)" proof (cases s)
            assume 1: "wf_row s (Suc m)" and 2: "s = []" 
            from 1 2 show ?thesis by (rule wf_help)
		next
            fix b listb 
			assume 5: "s = b # listb"
			from 2 4 have 6:"wf_row list m" by force
			from 3 5 have 7:"wf_row listb m" by force
			from 1 6 7 have 8:"wf_row (row_add list listb) m" by blast
			from 4 5 8 show "wf_row (row_add r s) (Suc m)"
				unfolding row_add_def by auto
		qed
	qed 
qed

lemma matrix_add_wf: "\<lbrakk> wf_matrix A m n; wf_matrix B m n \<rbrakk> \<Longrightarrow> wf_matrix (matrix_add A B) m n"
apply (induct m arbitrary: A B n)
apply (case_tac n)
apply clarify
apply (erule wf_matrix_null)+
apply clarify
unfolding matrix_add_def
apply (auto simp: row_add_wf)
apply (erule wf_matrix_null2)+
apply clarify
apply (auto simp: row_add_wf)
proof -
  fix m A B n
  assume 1: "\<And>A B n. \<lbrakk>wf_matrix A m n; wf_matrix B m n\<rbrakk> \<Longrightarrow> wf_matrix (map (\<lambda>(x, y). row_add x y) (zip A B)) m n"
     and 2: "wf_matrix A (Suc m) n"
     and 3: "wf_matrix B (Suc m) n"
  from 1 2 3 show "wf_matrix (map (\<lambda>(x, y). row_add x y) (zip A B)) (Suc m) n" apply (cases A) apply (cases B) apply force apply force apply (cases B) apply force 
  proof -
    fix a list aa lista
    assume 10: "\<And>A n B. \<lbrakk>wf_matrix A m n; wf_matrix B m n\<rbrakk> \<Longrightarrow> wf_matrix (map (\<lambda>(x, y). row_add x y) (zip A B)) m n"
       and 11: "wf_matrix A (Suc m) n"
       and 12: "wf_matrix B (Suc m) n"
       and 13: "A = a # list"
       and 14: "B = aa # lista"
    from 11 13 have 15: "wf_matrix list m n" by auto
    from 12 14 have 16: "wf_matrix lista m n" by auto
    from 10 15 16 have 17: "wf_matrix (map (\<lambda>(x, y). row_add x y) (zip list lista)) m n" by auto
    from 11 13 have 18: "wf_row a n" by auto
    from 12 14 have 19: "wf_row aa n" by auto
    from 17 13 14 18 19 show "wf_matrix (map (\<lambda>(x, y). row_add x y) (zip A B)) (Suc m) n" by (auto simp: row_add_wf)
  qed
qed

lemma add_equivalence: "\<lbrakk>matrix_equivalence mathA btoA m n; matrix_equivalence mathB btoB m n\<rbrakk>
       \<Longrightarrow> matrix_equivalence (math_add mathA mathB) (matrix_add btoA btoB) m n"
apply (erule matrix_equivalence.cases)
apply (auto)
apply (erule matrix_equivalence.cases)
apply (auto)
apply (erule matrix_add_wf, assumption)
apply (auto simp:add_access)
done


(*
------------
Matrix Multiplication
------------
*)


(* Well-Formedness *)

lemma fold_zip_wf: "\<lbrakk>wf_row v0 n; wf_matrix B m n; wf_row u m \<rbrakk> 
	\<Longrightarrow> wf_row (foldr (% (e,v) a . row_add (row_scale e v) a) (zip u B) v0) n"
apply (induct m arbitrary: v0 n B u e0)
apply (erule wf_matrix_null3)
apply simp
proof -
	fix m v0 n B u
    assume 1:"\<And>v0 n B u. \<lbrakk>wf_row v0 n; wf_matrix B m n; wf_row u m\<rbrakk> \<Longrightarrow> wf_row (foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip u B) v0) n"
	and 2:"wf_row v0 n" 
	and 3: "wf_matrix B (Suc m) n"
    and 4:"wf_row u (Suc m)"
	from 1 2 3 4 show "wf_row (foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip u B) v0) n"
	apply (cases B)	apply simp apply (cases u) apply simp
    proof -
		fix b0 subB u0 subu
		assume 5: "B = b0 # subB"
		and 6: "u = u0 # subu"
		from 3 5 have 7: "wf_matrix subB m n" by blast
		from 4 6 have 8: "wf_row subu m" by blast
		from 1 2 7 8 have 9: "wf_row (foldr (\<lambda>(e,v). row_add (row_scale e v)) (zip subu subB) v0) n" by blast
		from 3 5 have 10:"wf_row b0 n" by blast
		from 5 6 9 10 show "wf_row (foldr (\<lambda>(e,v). row_add (row_scale e v)) (zip u B) v0) n" by (auto simp: row_add_wf row_scale_wf)
	qed
qed


lemma vec_mat_wf: " \<lbrakk>wf_matrix B n p; wf_row a n; n > 0\<rbrakk> \<Longrightarrow> wf_row (vector_matrix_multiply a B) p"
apply (cases n)
apply auto
proof -
	fix nat
	assume 1:"wf_matrix B (Suc nat) p" 
	and 2:"wf_row a (Suc nat)"
	from 1 2 show "wf_row (vector_matrix_multiply a B) p" 
    apply (erule_tac A="B" and m="nat" and n="p" in wf_matrix_succ) apply (erule wf_row_succ) apply simp 
	apply (rule fold_zip_wf) by (auto simp: row_scale_wf)
qed

lemma matrix_mult_wf: "\<lbrakk> wf_matrix A m n; wf_matrix B n p; n > 0 \<rbrakk> \<Longrightarrow> wf_matrix (matrix_multiply A B) m p"
apply (induct m arbitrary: A B n p)
apply (erule wf_matrix_null3)
unfolding matrix_multiply_def
apply force
proof -
	fix m A B n p
    assume 1:"\<And>A B n p. \<lbrakk>wf_matrix A m n; wf_matrix B n p; 0 < n\<rbrakk> \<Longrightarrow> wf_matrix (map (\<lambda>u. vector_matrix_multiply u B) A) m p" 
	and 2:"wf_matrix A (Suc m) n" 
	and 3:"wf_matrix B n p"
	and ngtz:"n > 0"
	from 1 2 3 ngtz show "wf_matrix (map (\<lambda>u. vector_matrix_multiply u B) A) (Suc m) p" apply (cases A) apply (erule wf_matrix_succ) apply force 
	proof -
		fix a list
		assume 4: "A = a # list"  
		from 2 4 have 5:"wf_matrix list m n" by auto
		from 1 5 3 ngtz have 6:"wf_matrix (map (%u. vector_matrix_multiply u B) list) m p" by blast
		from 2 3 4 6 ngtz show ?thesis by (auto simp:vec_mat_wf)
	qed
qed

(* Value Equivalence *)

lemma shiftaccessmap: "\<lbrakk> a = x # suba; B = r # subB \<rbrakk> \<Longrightarrow>
	map (\<lambda> k. a ! k * B ! k ! j) [Suc 0..<Suc n] = map (\<lambda> k. suba ! k * subB ! k ! j) [0..<n]"
apply (induct n)
apply auto
done

lemma fold_sum_unroll: "\<lbrakk> a = x # suba; B = r # subB ; length a = (Suc n); length B = (Suc n)\<rbrakk> \<Longrightarrow>
	foldr op + (map (\<lambda>k. a ! k * B ! k ! j) [0..<Suc n]) (0::nat) = (x * (r ! j)) + (foldr op + (map (\<lambda>k. suba ! k * subB ! k ! j) [0..<n]) 0)"
proof - 
	assume 1:"a = x # suba" 
	and 2:"B = r # subB" 
	and 3:"length a = Suc n"
	and 4: "length B = Suc n"
	from 1 2 have 5:"foldr op + (map (\<lambda>k. a ! k * B ! k ! j) [0..<Suc n]) (0::nat) = 
		(x * (r ! j)) + (foldr op + (map (\<lambda>k. a ! k * B ! k ! j) [Suc 0..<Suc n]) (0::nat))" by (simp add: upt_conv_Cons foldr_map not_le)
	from 1 2 also have 6:"\<dots> = (x  * (r ! j)) + (foldr op + (map (\<lambda>k. suba ! k * subB ! k ! j) [0..<n]) 0)"
		apply (rule_tac s="(map (\<lambda>k. a ! k * B ! k ! j) [Suc 0..<Suc n])" and t="(map (\<lambda>k. suba ! k * subB ! k ! j) [0..<n])" in subst)
		apply (erule shiftaccessmap) apply assumption by blast
	from 5 6 show ?thesis by simp
qed

(* Key lemma: vector_matrix multiply gives correct value *)
lemma row_mat_access: "\<lbrakk> wf_matrix B n p; wf_row a n; j < p; n > 0 \<rbrakk> \<Longrightarrow>
	vector_matrix_multiply a B ! j = foldr op + (map (\<lambda>k. a ! k * B ! k ! j) [0..<n]) 0"
(* First we need to case tac a & B.  This corresponds to the cases in the def of vector_matrix_multiply.
Neither can be [] because n > 0.  We case_tac n (the size) for the same reason, to get to n = Suc size
*)
apply (cases B) apply clarify apply (erule wf_matrix_empty) apply simp (* proves B != [] *)
apply (cases a) apply clarify apply (erule wf_row_empty) apply simp (* proves a != [] *)
apply (cases n) apply clarify (* proves n != 0 *)
(* Now we know a & B & n all have sub-parts.  We fix them, then we can induct size (the subpart of n *)
thm subst
apply (rule_tac s="Suc nat" and t="n" in subst) apply clarify
thm ssubst
apply (rule_tac t="foldr op + (map (\<lambda>k. a ! k * B ! k ! j) [0..<Suc nat]) 0" and s="(aaa * (aa ! j)) + (foldr op + (map (\<lambda>k. lista ! k * list ! k ! j) [0..<nat]) 0)" in ssubst)
apply (rule fold_sum_unroll) apply simp apply simp apply simp apply (rule wf_row_implies_len) apply (erule wf_row2) apply assumption
thm wf_matrix_implies_len
apply (rule_tac n="p" in wf_matrix_implies_len) apply clarify
proof clarify
	fix B0 B' a0 a' size
	assume 1: "j < p"
	and 2: "wf_row B0 p"
	and 3: "wf_matrix B' size p"
	and 4: "wf_row a' size"
	from 1 2 3 4 have proof_form:"foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j =
        a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<size]) 0"
	proof (induct size arbitrary: j p B0 B' a0 a')
		(* Now we enter the meat of the proof: inducting on size. *)
		(* First the base case: size = 0.  In this case a' and B' will both be 0 *)
		fix j p B0 B' a0 a'
		assume 1: "j < p"
		and 2: "wf_row B0 p"
		and 3: "wf_matrix B' 0 p"
		and 4: "wf_row a' 0"
		from 3 have 5:"B' = []" apply (rule wf_matrix_null3) by assumption
		from 4 have 6:"a' = []" apply (rule wf_row_null) by assumption
		from 1 2 5 6 show "foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j = 
			a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<0]) 0" by (simp add: scale_row_access)
	next
		fix size j p B0 B' a0 a'
                assume 1: "\<And>j p B0 B' a0 a'. \<lbrakk>j < p; wf_row B0 p; wf_matrix B' size p; wf_row a' size\<rbrakk>
                            \<Longrightarrow> foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j = 
				a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<size]) 0"
		and 2: "j < p" 
		and 3: "wf_row B0 p" 
		and 4: "wf_matrix B' (Suc size) p" 
		and 5: "wf_row a' (Suc size)"
		(* Now we need to again use cases on B' and a' to get at their sub-parts *)
		from 1 2 3 4 5 show "foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j = 
			a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<(Suc size)]) 0" 
		apply (cases B') apply (erule wf_matrix_succ) apply simp
		apply (cases a') apply (erule wf_row_succ) apply simp
		proof -
		fix B1 B'' a1 a''
		assume 6:"B' = B1 # B'' "
		assume 7:"a' = a1 # a'' "
		(* Here we gather up all the pieces we need *)
		from 4 6 have 8: "wf_matrix B'' size p" apply simp apply (erule wf_matrix2) by assumption
		from 5 7 have 9: "wf_row a'' size" apply simp apply (erule wf_row2) by assumption
		from 4 6 have 12: "wf_row B1 p" apply simp apply (erule wf_matrix2) by assumption
		from 12 have 13: "wf_row (row_scale a1 B1) p" by (rule row_scale_wf)
		from 3 have 14: "wf_row (row_scale a0 B0) p" by (rule row_scale_wf)
		from 8 9 14 have 15: "wf_row (foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a'' B'') (row_scale a0 B0)) p" by (simp add:fold_zip_wf)
		from 4 have 18:"length B' = Suc size" by (rule wf_matrix_implies_len)
		from 5 have 19:"length a' = Suc size" by (rule wf_row_implies_len)
		from 2 12 have 22:"(row_scale a1 B1) ! j = a1 * B1 ! j" by (simp add: scale_row_access)
        from 7 6 19 18 have 20: "foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<Suc size]) (0::nat) = 
			(a1 * (B1 ! j)) + (foldr op + (map (\<lambda>k. a'' ! k * B'' ! k ! j) [0..<size]) 0)" by (rule fold_sum_unroll)
		(* And finally step through the chain of equations to show equality *)
		from 6 7 have 11:"foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j =
			(row_add (row_scale a1 B1) (foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a'' B'') (row_scale a0 B0))) ! j" by auto
		also from 2 13 15 have 
			"\<dots> = (row_scale a1 B1) ! j + (foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a'' B'') (row_scale a0 B0)) ! j" by (simp add: add_row_access)
		also from 1 2 3 8 9 have
			 "\<dots> = (row_scale a1 B1) ! j + (a0 * B0 ! j + foldr op + (map (\<lambda>k. a'' ! k * B'' ! k ! j) [0..<size]) 0)" by auto
		also from 22 have 
			"\<dots> = a1 * B1 ! j + (a0 * B0 ! j + foldr op + (map (\<lambda>k. a'' ! k * B'' ! k ! j) [0..<size]) 0)" by simp
		also from 20 have 
			"\<dots> = a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<Suc size]) (0::nat)" by simp
		finally show "foldr (\<lambda>(e, v). row_add (row_scale e v)) (zip a' B') (row_scale a0 B0) ! j = 
			a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<Suc size]) 0" by simp
	qed qed
	from proof_form show "vector_matrix_multiply (a0 # a') (B0 # B') ! j = 
		a0 * B0 ! j + foldr op + (map (\<lambda>k. a' ! k * B' ! k ! j) [0..<size]) 0" by simp
qed


lemma matrix_mult_access: "\<lbrakk> wf_matrix A m n; wf_matrix B n p; i < m; j < p; n > 0 \<rbrakk> \<Longrightarrow> 
	matrix_access (matrix_multiply A B) i j = foldr plus (map (% k . (matrix_access A i k) * (matrix_access B k j)) [0..<n]) 0"
apply (induct m arbitrary: A n B p i j)
apply auto
proof - 
	fix m A n B p i j
       assume 1: "\<And>A n B p i j.
           \<lbrakk>wf_matrix A m n; wf_matrix B n p; i < m; j < p; 0 < n\<rbrakk>
           \<Longrightarrow> matrix_access (matrix_multiply A B) i j = foldr op + (map (\<lambda>k. matrix_access A i k * matrix_access B k j) [0..<n]) 0"
        and 2: "wf_matrix A (Suc m) n" 
	and 3: "wf_matrix B n p" 
	and 4: "i < Suc m" 
	and 5: "j < p"
	and 6: "0 < n"
    from 1 2 3 4 5 show "matrix_access (matrix_multiply A B) i j = foldr op + (map (\<lambda>k. matrix_access A i k * matrix_access B k j) [0..<n]) 0"
	apply (cases A) apply (erule wf_matrix_succ) apply force
	proof -
		fix a list
		assume 7: "A = a # list"
		from 2 7 have 8: "wf_matrix list m n" by auto
		from 1 3 4 5 6 8 show ?thesis
		proof (cases i)
			assume 9:"i=0"	
			from 1 2 3 4 5 6 7 8 9 show ?thesis unfolding matrix_access_def unfolding matrix_multiply_def apply simp apply (erule row_mat_access) by auto
		next
			fix nat
			assume 9: "i = Suc nat"
			from 7 9 have 12: "matrix_access (matrix_multiply A B) i j = matrix_access (matrix_multiply list B) nat j" 
				unfolding matrix_access_def unfolding matrix_multiply_def by auto
			from 9 4 have 10:"nat < m" by simp
			from 1 8 3 10 5 6 have 
				11:"matrix_access (matrix_multiply list B) nat j = foldr op + (map (\<lambda>k. matrix_access list nat k * matrix_access B k j) [0..<n]) 0" by simp
			from 9 7 have 13: "foldr op + (map (\<lambda>k. matrix_access A i k * matrix_access B k j) [0..<n]) 0 = 
				foldr op + (map (\<lambda>k. matrix_access list nat k * matrix_access B k j) [0..<n]) 0" unfolding matrix_access_def by auto
			from 12 11 13 show ?thesis by auto
		qed
	qed
qed

lemma mult_equivalence: "\<lbrakk> matrix_equivalence mathA btoA m n; matrix_equivalence mathB btoB n p; n > 0\<rbrakk>
       \<Longrightarrow> matrix_equivalence (math_mult mathA mathB) (matrix_multiply btoA btoB) m p"
apply (erule matrix_equivalence.cases)
apply auto
apply (erule matrix_equivalence.cases)
apply auto
apply (erule matrix_mult_wf, assumption)
apply (auto simp:matrix_mult_access)
apply (rule foldr_cong)
apply auto
done

lemma wf_expression_greater_than_zero: "wf_expression e m n \<Longrightarrow> n > 0"
apply (induct rule: wf_expression.induct)
apply auto
done

(* Main Theorem *)

theorem compiler_equivalence: "wf_expression e m n \<Longrightarrow> matrix_equivalence (math_eval e) (compute e) m n"
apply (induct rule: wf_expression.induct)
apply auto
apply (erule add_equivalence, assumption)
apply (erule mult_equivalence, assumption)
apply (erule wf_expression_greater_than_zero)
apply (erule scale_equivalence)
done













(* Example *)
(* 
1 2
3 4
5 6
*
1 0
4 4
=
9 8
19 16
29 24
*)

definition example :: expr where
  "example = Mult (Matrix [[1, 2], [3, 4], [5, 6]] 3 2) (Matrix [[1, 0], [4, 4]] 2 2)"

value "matrix_access (compute example) 1 0"
value "math_matrix_fun (math_eval example) (1, 0)" 

value "(matrix_access (compute example) 1 0) = (math_matrix_fun (math_eval example) (1, 0))"

end
